#include "luaclass.h"
#include "cpgplot.h"
#include "fitsio.h"
#include <iostream>
#include <fstream>
#include <cmath>
#include <sys/types.h>
#include <sys/stat.h>
#include <unistd.h>
#include <cstring>

#include "eclipse_types.h"

//Lua lua("");

std::vector<Image *> images;
int Fits2Float (const char *, float **array, int *arx, int *ary);

extern float pixelscale_x;
extern float pixelscale_y;
extern float xref;
extern float yref;
extern int sizex;
extern int sizey;
extern int xmin;
extern int xmax;
extern int ymin;
extern int ymax;

extern float zoom;
extern float* data;

extern bool verbose;

extern std::vector<int> cont_colors;
extern std::vector<int> cont_styles;
extern std::vector<float> cont_values;

float ContourSampling = 1.0;
float ContourSmoothing = 0.0;
float ContourStart = 2.5;
bool ContourSaving = true;

void pgellipse_c( double Cpx, 
                  double Cpy, 
                  double Major, 
                  double Minor, 
                  double Pa,
                  double Startangle,
                  double Endangle,
                  double Delta );
void pgcone_c( double Cpx, 
	       double Cpy, 
	       double R, 
	       double Startangle,
	       double Endangle,
	       double Delta );
void DefaultContours (float *, int, int, std::vector<float> *);
float *  LinearSampling (float *array, int nx, int ny, double*trans, int *nnx, int *nny);
float *Smooth (float *array, int nx, int ny, double sigma);
float *FastSmooth (float *array, int nx, int ny, double sigma);
void SaveFits (float *array, int nx, int ny, const char*name);

class CallNode {
public:
  const char * name;
  lua_CFunction func;
  int status;
};

CallNode calls[] = {
  {"Drawing primitives", NULL,0},
  {"DrawFilledRectangle", Lua::FilledRectangle,0},
  {"DrawRectangle", Lua::Rectangle, 0},
  {"Rectangle",     Lua::Rectangle, 1},
  {"DrawEllipse",   Lua::Ellipse,   0},
  {"DrawLine",      Lua::DrawLine,  0},
  {"DrawArrow",     Lua::DrawArrow, 0},
  {"DrawText",      Lua::Text,      0},
  {"DrawCone",      Lua::DrawCone,  0},
  {"DrawCross",     Lua::DrawCross, 0},
  {"Text",          Lua::Text,      1},
  
  {"Other drawing functions", NULL, 0},
  {"Contour",  Lua::Contour,        0},
  {"DrawSlit",  Lua::DrawSlit,      0},
  {"DrawRegion",  Lua::DrawRegion,  0},
  {"DrawScaleBar", Lua::DrawScaleBar, 0},
  {"DrawScalebar", Lua::DrawScaleBar, 1},

  {"Output variables", NULL,        0},
  {"SetTextSize",  Lua::TextSize,   0},
  {"TextSize",     Lua::TextSize,   1},
  {"SetLineWidth", Lua::LineWidth,  0},
  {"LineWidth",    Lua::LineWidth,  1},
  {"SetColor",     Lua::SetColor,   0},
  {"SetColour",    Lua::SetColor,   1},
  {"SetLineStyle", Lua::LineStyle,  0},
  {"SetFillStyle", Lua::SetFillStyle,0},
  {"SetContourSampling", Lua::SetContourSampling,0},
  {"SetContourSmoothing", Lua::SetContourSmoothing,0},
  {"SetContourSaving", Lua::SetContourSaving,0},
    
  {"Data related I/O primitives", NULL,0},
  {"LoadImage",  Lua::LoadImage,       0},
  {"SaveFits",   Lua::SaveFits,        0},
  {"FreeImage",  Lua::FreeImage,       0},
  {"PxReplace",  Lua::PxReplace,       0},
  {"PxArith",  Lua::PxArith,           0},
  {"ClipImage", Lua::ClipImage,        0},

  {"Misc Functions", NULL,             0},
  {"Input",  Lua::Input,               0},
  {"ListCommands", Lua::ListCommands,  0},

  {"API to pgplot calls", NULL,             0},
  {"ViewPort", Lua::ViewPort, 0},
  {"Viewport", Lua::ViewPort, 1},
  {"Limits", Lua::Limits, 0},
  {"Box", Lua::Box, 0},
  {"Plot", Lua::Plot, 0},
  {NULL, NULL,0}
};


bool Lua::ExecString (const char *str)
{
  //printf ("%s\n",str);
  int error = luaL_loadbuffer(L, str, strlen(str), str) || lua_pcall(L, 0, 0, 0);
  return (error);
}

Lua::Lua (const char *file)
{
  L = lua_open();

  /* load Lua base libraries */
  //lua_baselibopen(L);
  luaL_openlibs(L);
  
  /* register functions to be used in lua scripts */
  unsigned int i=0;
  while (calls[i].name != NULL) {
    if (calls[i].func != NULL)
      lua_register (L, calls[i].name, calls[i].func);
    i++;
  }


//   lua_register(L, "SubSample",  SubSample);
//   lua_register(L, "SmoothGaussian", SmoothGaussian);
//   lua_register(L, "ClipImage", ClipImage);

  // lua_register(L, "Region", Region);

  /* run the script */
  if (file != NULL)
    // changed
    luaL_dofile(L, file);
}

Lua::~Lua ()
{
  /* cleanup Lua */
  lua_close(L); 

  for (unsigned int i=0; i<images.size(); i++)
    delete images[i];
}

int Lua::DrawRegion (lua_State *L)
{
  CheckArgsNumericAlpha ch(L, 0, 1);
  if (!ch.Success()) {
    printf ("USAGE: DrawRegion (\"file\");\n");
    return (MethodError (L, "DrawRegion"));
  }
  int scale=1;

  if (zoom>-1)
    scale = (int)zoom;

  const char *chs = lua_tostring (L,1);

  // printf ("DrawRegion: creating image with size of %dx%d\n",sizex,sizey);
  // printf ("DrawRegion: xref=%f, yref=%f\n", xref, yref);
  float *data = new float[sizex*sizey];

  for (int i=0; i<sizex*sizey; i++) data[i]=0.0;

  FILE *file=fopen (chs, "r");
  while (!feof(file)) {
    int x,y;
    if (fscanf (file, " %d %d", &x, &y)!=2)
      break;
    // printf ("%d %d\n", x, y);
    if (x>=xmin & x<=xmax & y>=ymin & y<=ymax)
      for (int i=0; i<scale; i++)
	for (int j=0; j<scale; j++)
	  data[scale*(x-xmin)+i+(scale*(y-ymin)+j)*sizex] = 1.0;
  }

  float cont[2] = {0.5, 1.5};

  //  float tr[6] = {xref, -pixelscale_x, 0.0, -yref, 0.0, pixelscale_y};
  float tr[6] = {xref, -pixelscale_x, 0.0, -yref, 0.0, pixelscale_y};

  cpgcont (data, sizex, sizey, 
	   1, sizex, 1, sizey, 
	   &cont[0], -2, tr);

  delete[] data;

  return (Success(L));


}

bool Lua::ExecFile (const char *file)
{
    // changed
  return (!luaL_dofile(L, file));
}

int Lua::Input (lua_State *L)
{
  CheckArgsNumericAlpha ch(L, 0, 1);
  if (!ch.Success()) {
    printf ("USAGE: Input (\"file\");");
    return (MethodError (L, "Input"));
  }

    // changed
  return (luaL_dofile(L, ch.text[0]));
}

int Lua::Rectangle (lua_State *L)
{
  CheckArgsNumeric ch(L, 4);
  if (!ch.Success()) {
    printf ("USAGE: DrawRectangle (x1, y2, x2, y2);\n");
    return (MethodError (L, "Rectangle"));
  }

  int n = 5;
  float x[5] = {ch.val[0], ch.val[2], ch.val[2], ch.val[0], ch.val[0]};
  float y[5] = {ch.val[1], ch.val[1], ch.val[3], ch.val[3], ch.val[1]};
  cpgline (n, x, y);
  return (Success(L));
}

int Lua::FilledRectangle (lua_State *L)
{
  CheckArgsNumeric ch(L, 4);
  if (!ch.Success()) {
    printf ("USAGE: DrawFilledRectangle (x1, y2, x2, y2);\n");
    return (MethodError (L, "DrawFilledRectangle"));
  }

  cpgrect (ch.val[0], ch.val[2],ch.val[1], ch.val[3] );
  return (Success(L));
}

int Lua::ListCommands (lua_State *L)
{
  int i=0;
  int j;
  while (calls[i].name != NULL) {
    j=0;
    bool depr = false;
    if (calls[i].func == NULL) 
      printf("\n%s:\n", calls[i].name);
    else {
      if (calls[i].status == 2)
	printf ("  %s (deprecated)\n", calls[i].name);
      else if (calls[i].status == 0)
	printf ("  %s\n", calls[i].name);
    }

//       while (j<i) {
// 	if (calls[i].func == calls[j].func) {
// 	  depr = true;
// 	  break;
// 	}
// 	j++;
//       }
//       if (depr)
// 	printf ("  %s (deprecated)\n", calls[i].name);
//       else
// 	printf ("  %s\n", calls[i].name);
      
//     }
    i++;
  }
  printf ("\n");
  return (Success(L));
}

int Lua::DrawLine (lua_State *L)
{
  CheckArgsNumeric ch(L, 4);
  if (!ch.Success()) {
    printf ("USAGE: DrawLine(x1, y1, x2, y2);\n");
    return (MethodError (L, "DrawLine"));
  }

  int n = 2;
  float x[2] = {ch.val[0], ch.val[2]};
  float y[2] = {ch.val[1], ch.val[3]};
  cpgline (n, x, y);
  return (Success(L));
}

int Lua::DrawArrow (lua_State *L)
{
  if (lua_gettop(L) == 4)
    return (DrawArrow4 (L));

  if (lua_gettop(L) == 7)
    return (DrawArrow7 (L));

  printf ("USAGE: DrawArrow (x1, y1, x2, y2);\n");
  printf ("       DrawArrow (x1, y1, x2, y2, fs, a, b);\n");
  printf ("       where fs: 1 (filled), 2(outline) (default 1)\n");
  printf ("             a: angle (default 45)\n");
  printf ("             a: barb (default 0.3)\n");
  return (MethodError (L, "DrawLine Sel."));
}


int Lua::DrawArrow4 (lua_State *L)
{
  CheckArgsNumeric ch(L, 4);
  if (!ch.Success()) {
    printf ("USAGE: DrawArrow (x1, y1, x2, y2);\n");
    return (MethodError (L, "DrawArrow4"));
  }

  cpgarro (ch.val[0], ch.val[1], ch.val[2], ch.val[3]);
  return (Success(L));
}

int Lua::DrawArrow7 (lua_State *L)
{
  CheckArgsNumeric ch(L, 7);
  if (!ch.Success()) {
    printf ("USAGE: DrawArrow (x1, y1, x2, y2, fs, a, b);\n");
    printf ("       where fs: 1 (filled), 2(outline) (default 1)\n");
    printf ("             a: angle (default 45)\n");
    printf ("             a: barb (default 0.3)\n");
    return (MethodError (L, "DrawArrow7"));
  }

  if (ch.val[4]>1.5)
    cpgsah (2, ch.val[5], ch.val[6]);
  else
    cpgsah (1, ch.val[5], ch.val[6]);
   
  cpgarro (ch.val[0], ch.val[1], ch.val[2], ch.val[3]);

  return (Success(L));
}

int Lua::Ellipse (lua_State *L)
{
  float cx, cy, minor, major, pa, pa_start,pa_end;
  if (lua_gettop(L)==5) {
    CheckArgsNumeric ch(L, 5);
    if (!ch.Success()) {
      printf ("USAGE: x y minor major pa\n");
      return (MethodError (L, "DrawEllipse 1"));
    }

    cx    = ch.val[0];
    cy    = ch.val[1];
    minor = ch.val[2];
    major = ch.val[3];
    pa    = ch.val[4];
    pa_start = 0.0;
    pa_end = 360.0;
  }
  else if (lua_gettop(L)==7) {
    CheckArgsNumeric ch(L, 7);
    if (!ch.Success()) {
      printf ("USAGE: x y minor major pa pa_start pa_end\n");
      return (MethodError (L, "DrawEllipse 2"));
    }

    cx    = ch.val[0];
    cy    = ch.val[1];
    minor = ch.val[2];
    major = ch.val[3];
    pa    = ch.val[4];
    pa_start = ch.val[5];
    pa_end = ch.val[6];
  }
  else {
    printf ("USAGE: x y minor major pa\n");
    printf ("       x y minor major pa pa_start pa_end\n");
    return (MethodError (L, "DrawEllipse 3"));
  }

  pgellipse_c (cx, cy, minor, major, pa, pa_start, pa_end, 1.0);

  return (Success(L));
}

int Lua::DrawScaleBar (lua_State *L)
{
  CheckArgsNumeric ch(L, 3);
  if (!ch.Success()) {
    printf ("USAGE: x y length\n");
    return (MethodError (L, "DrawScaleBar"));
  }

  cpgmove (ch.val[0]-0.5*ch.val[2], ch.val[1]-0.07*ch.val[2]);
  cpgdraw (ch.val[0]-0.5*ch.val[2], ch.val[1]+0.07*ch.val[2]);

  cpgmove (ch.val[0]+0.5*ch.val[2], ch.val[1]-0.07*ch.val[2]);
  cpgdraw (ch.val[0]+0.5*ch.val[2], ch.val[1]+0.07*ch.val[2]);
  
  cpgmove (ch.val[0]-0.5*ch.val[2], ch.val[1]);
  cpgdraw (ch.val[0]+0.5*ch.val[2], ch.val[1]);

  return (Success(L));
}

int Lua::DrawCross (lua_State *L)
{
  CheckArgsNumeric ch(L, 3);
  if (!ch.Success()) {
    printf ("USAGE: x y size\n");
    return (MethodError (L, "DrawCross"));
  }

  cpgmove (ch.val[0]-0.5*ch.val[2], ch.val[1]-0.5*ch.val[2]);
  cpgdraw (ch.val[0]+0.5*ch.val[2], ch.val[1]+0.5*ch.val[2]);

  cpgmove (ch.val[0]-0.5*ch.val[2], ch.val[1]+0.5*ch.val[2]);
  cpgdraw (ch.val[0]+0.5*ch.val[2], ch.val[1]-0.5*ch.val[2]);
  
  return (Success(L));
}


int Lua::DrawCone (lua_State *L)
{
  CheckArgsNumeric ch(L, 5);
  if (!ch.Success()) {
    printf ("USAGE: x y r pa width\n");
    return (MethodError (L, "DrawCone"));
  }

  float cx    = ch.val[0];
  float cy    = ch.val[1];
  float r = ch.val[2];
  float pa = ch.val[3];
  float width    = ch.val[4];

  pgcone_c (cx, cy, r, pa-width, pa+width, 1.0);

  return (Success(L));

}

int Lua::DrawSlit (lua_State *L)
{
  CheckArgsNumeric ch(L, 3);
  if (!ch.Success()) {
    printf ("USAGE: DrawSlit (pa, half_length, half_width);\n");
    return (MethodError (L, "DrawSlit"));
  }

  int n = 5;
  
  float pa = ch.val[0];
  float hlength = ch.val[1]/2.0;
  float hwidth = ch.val[2]/2.0;

  float cpa = std::cos (pa*0.01745329251);
  float spa = std::sin (pa*0.01745329251);
  
  float x[5];
  float y[5];

  x[0] = -( hwidth*cpa-hlength*spa);
  x[1] = -( hwidth*cpa+hlength*spa);
  x[2] = -(-hwidth*cpa+hlength*spa);
  x[3] = -(-hwidth*cpa-hlength*spa);
  x[4] = -( hwidth*cpa-hlength*spa);

  y[0] = hwidth*spa+hlength*cpa;
  y[1] = hwidth*spa-hlength*cpa;
  y[2] =-hwidth*spa-hlength*cpa;
  y[3] =-hwidth*spa+hlength*cpa;
  y[4] = hwidth*spa+hlength*cpa;

  cpgline (n, x, y);
  return (Success(L));
}

int Lua::TextSize (lua_State *L)
{
  CheckArgsNumeric ch(L, 1);
  if (!ch.Success()) {
    printf ("USAGE: TextSize (2);\n");
    return (MethodError (L, "TextSize"));
  }

  cpgsch (ch.val[0]);
  return (Success(L));
}

int Lua::LineWidth (lua_State *L) 
{
  CheckArgsNumeric ch(L, 1);
  if (!ch.Success()) {
    printf ("USAGE: LineWidth (2);\n");  
    return (MethodError (L, "LineWidth"));
  }

  cpgslw ((int)(ch.val[0]+0.5));
  return (Success(L));
}

int Lua::Text (lua_State *L)
{
  CheckArgsNumericAlpha ch(L, 4, 1);
  if (!ch.Success()) {
    printf ("USAGE: DrawText (x, y, angle, justi, \"text\");\n");
    printf ("       justification:\n");
    printf ("         0.0: left, 0.5 centered, right 1.0\n");
    return (MethodError (L, "Text"));
  }
  if (verbose)
    printf ("Luaclass::Text \"%s\"\n", ch.text[0]);

  cpgptxt (ch.val[0], ch.val[1], ch.val[2], ch.val[3], ch.text[0]);
  return (Success(L));
}

int Lua::SetColor (lua_State *L)
{
  CheckArgsNumeric ch(L, 1);
  if (!ch.Success())
    return (MethodError (L, "SetColor"));

  cpgsci ((int)(ch.val[0]+0.5));
  return (Success(L));
}

int Lua::SetContourSampling (lua_State *L)
{
  CheckArgsNumeric ch(L, 1);
  if (!ch.Success())
    return (MethodError (L, "SetContourSampling"));

  ContourSampling = ch.val[0];
  return (Success(L));
}


int Lua::SetContourSmoothing (lua_State *L)
{
  CheckArgsNumeric ch(L, 1);
  if (!ch.Success())
    return (MethodError (L, "SetContourSmoothing"));

  ContourSmoothing = ch.val[0];
  return (Success(L));
}

int Lua::SetContourSaving (lua_State *L)
{
  CheckArgsAlphaNumeric ch(L, 1, 0);
  if (!ch.Success())
    return (MethodError (L, "SetContourSaving"));

  if (strcmp(ch.text[0], "yes")==0)
    ContourSaving = true;
  else
    ContourSaving = false;
  return (Success(L));
}

int Lua::LineStyle(lua_State *L)
{
  CheckArgsNumeric ch(L, 1);
  if (!ch.Success()) {
    printf ("USAGE: SetLineStyle (l);\n");
    printf ("       1: full line\n       2: dashed\n       3: dot-dash-dot-dash\n");
    printf ("       4: dotted\n       5: dash-dot-dot-dot\n");
    return (MethodError (L, "SetLineStyle"));
  }

  if (ch.val[0]<0.5 || ch.val[0]>5.5) {
    printf ("USAGE: SetLineStyle (l);\n");
    printf ("       1: full line\n       2: dashed\n       3: dot-dash-dot-dash\n");
    printf ("       4: dotted\n       5: dash-dot-dot-dot\n");
    return (MethodError (L, "SetLineStyle"));
  }
  
  cpgsls ((int)(ch.val[0]+0.5));
  return (Success(L));
}

int Lua::SetFillStyle(lua_State *L)
{
  CheckArgsNumeric ch(L, 1);
  if (!ch.Success()) {
    printf ("USAGE: SetFillStyle (fs);\n");
    printf ("  FS = 1 => solid");
    printf ("  FS = 2 => outline");
    printf ("  FS = 3 => hatched");
    printf ("  FS = 4 => cross-hatched");
    return (MethodError (L, "SetFillStyle"));
  }

  if (ch.val[0]<0.5 || ch.val[0]>4.5) {
    printf ("USAGE: SetFillStyle (fs);\n");
    printf ("  FS = 1 => solid");
    printf ("  FS = 2 => outline");
    printf ("  FS = 3 => hatched");
    printf ("  FS = 4 => cross-hatched");
    return (MethodError (L, "SetFillStyle"));
  }
  
  cpgsfs ((int)(ch.val[0]+0.5));
  return (Success(L));
}

int Lua::ClipImage (lua_State *L)
{
  if (lua_gettop (L) == 2) {
    if ( !(lua_isnumber (L, 1) &&lua_isnumber (L, 2)))
      return (MethodError (L, "ClipImage1"));
    
    float lv = lua_tonumber (L, 1);
    float uv = lua_tonumber (L, 2);

    for (int i=0; i<sizex; i++)
      for (int j=0; j<sizey; j++) {
	if (data[i+j*sizex]<lv && data[i+j*sizex]!=-1e30)
	  data[i+j*sizex] = lv;
	if (data[i+j*sizex]>uv)
	  data[i+j*sizex] = uv;
      }
    return (Success(L));
  }
  if (lua_gettop (L) == 3) {
    CheckArgsAlphaNumeric ch(L, 1, 2);
    if (!ch.Success())
      return (MethodError (L, "ClipImage2"));

    int index = FindImage (ch.text[0]);
    if (index <0) {
      fprintf (stderr, "ERROR: Cannot find image %s\n", ch.text[0]);
      lua_pushnumber (L, 1);
      return (1);
    }

    for (unsigned int i=0; i<images[index]->sizex; i++)
      for (unsigned int j=0; j<images[index]->sizey; j++) {
	if (images[index]->data[i+j*sizex]<ch.val[0] 
	    && images[index]->data[i+j*sizex]!=-1e30)
	  images[index]->data[i+j*sizex] = ch.val[0];
	if (images[index]->data[i+j*sizex]>ch.val[1])
	  images[index]->data[i+j*sizex] = ch.val[1];
      }
    return (Success(L));
  }
	
  return (MethodError(L, "ClipImage3"));
}

int Lua::SaveFits (lua_State *L)
{
  fitsfile *fptr;
  int status = 0;
  CheckArgsAlphaNumeric ch(L, 2, 0);
  if (!ch.Success())
    return (MethodError (L, "SaveFits"));

  int index = FindImage (ch.text[0]);
  if (index <0) {
    std::cout << "ERROR: Cannot find image " << ch.text[0] 
	      << std::endl;
    lua_pushnumber (L, 1);
    return (1);
  }

  long naxes[2] = {images[index]->sizex, images[index]->sizey};
  long fpixel[2] = {1, 1};
  fits_create_file  (&fptr, ch.text[1], &status);
  if (status) {
    std::cout << "ERROR: Cannot create file " << ch.text[1]
	      << " - it probably exists already" 
	      << std::endl;
    lua_pushnumber (L, 1);
    return (1);
  }

  fits_create_img(fptr, FLOAT_IMG, 2, naxes, &status);

  fits_write_pix(fptr, TFLOAT, fpixel, 
		 images[index]->sizex*images[index]->sizey,
		 images[index]->data, &status);
  fits_close_file (fptr, &status);

  return (Success(L));
}

int Lua::Contour (lua_State *L)
{
  if (verbose)
    printf ("INFO: Contour\n");

  if (lua_gettop(L)!=2 && lua_gettop(L)!=1 && lua_gettop(L)!=5) {
    printf ("Usage: Contour (\"image.fits\", cx, cy, pxsize, {3.0, 3.0});\n");
    printf ("Usage: Contour (\"image\", {3.0, 3.0});\n");
    printf ("       Contour ({3.0, 3.0});\n");
    return (MethodError (L, "Contour"));
  }

  int curstyle, curcolor;
  cpgqls (&curstyle);
  cpgqci (&curcolor);

  int tindex = lua_gettop(L); // The pos. of table in stack
  std::vector<float> fff;
  lua_pushnil(L);  /* first key */
  while (lua_next(L, tindex) != 0) {
    /* `key' is at index -2 and `value' at index -1 */
    float temp = lua_tonumber(L, -1);
    fff.push_back (temp);
    if (ContourSaving) {
      cont_colors.push_back (curcolor);
      cont_styles.push_back (curstyle);
      cont_values.push_back (temp);
    }
    lua_pop(L, 1);  /* removes `value'; keeps `key' for next iteration */ 
  }

  if (fff.size()==0) {
    if (verbose) {
      printf ("Cannot get contour levels, will estimate in DefaultContours\n");
      printf ("This will likely OVERESTIMATE sigma!\n");
    }
  }

  double lineartrans[6];
  lineartrans[0] = ContourSampling ;
  lineartrans[1] = 0.00 ;
  lineartrans[2] = 0.00 ;
  lineartrans[3] = 0.00 ;
  lineartrans[4] = ContourSampling ;
  lineartrans[5] = 0.00 ;

  int _sizex, _sizey;
  float *_data;
  float _xref, _yref, _pixelscale_x, _pixelscale_y;
  if (lua_gettop(L) == 1) {
    _sizex = sizex;
    _sizey = sizey;
    _data = data;
    _xref = xref;
    _yref = yref;
    _pixelscale_x = pixelscale_x;
    _pixelscale_y = pixelscale_y;
  }
  else if (lua_gettop(L) == 2) {
    int index = FindImage (lua_tostring (L,1));
    if (index <0) {
      const char *ch = lua_tostring (L,1);
      std::cout << "ERROR: Cannot find image " << ch << std::endl;
      lua_pushnumber (L, 1);
      return (1);
    }
    float scale = (float)sizex/(float)images[index]->sizex;
    _sizex = images[index]->sizex;
    _sizey = images[index]->sizey;
    _data = images[index]->data;
    _xref = xref;
    _yref = yref;
    _pixelscale_x = pixelscale_x*scale;
    _pixelscale_y = pixelscale_y*scale;
  }
  else if (lua_gettop(L) == 5) {
    if (verbose)
      printf ("Contour5: %d %d\n", _sizex, _sizey);

    Fits2Float (lua_tostring (L,1), &_data, &_sizex, &_sizey);
	     
    _pixelscale_x = lua_tonumber (L, 4);
    _pixelscale_y = lua_tonumber (L, 4);
    _xref = (float)(lua_tonumber (L, 2)*_pixelscale_x);
    _yref = (float)(lua_tonumber (L, 3)+1.0)*_pixelscale_y;
    printf ("%f %f\n", _xref, _xref);
  }
  else {
    printf ("Usage: Contour (\"image.fits\", cx, cy, pxsize, {3.0, 3.0});\n");
    printf ("Usage: Contour (\"image\", {3.0, 3.0});\n");
    printf ("       Contour ({3.0, 3.0});\n");
    return (MethodError (L, "Contour"));
  }

  //::SaveFits (_data, _sizex, _sizey, "contour.fits");

//   if (lua_gettop(L) == 1) {
  if (ContourSampling != 1.0 && ContourSmoothing != 0.0) {
    double sigma = ContourSmoothing;
    
    if (verbose)
      printf ("Contour 0A: Smoothing with sigma=%.2f pixels\n", sigma);
    float *newdata =  FastSmooth (_data, _sizex, _sizey, sigma);
    //::SaveFits (_data, _sizex, _sizey, "contour0a.fits");
    if (fff.size()==0)
      DefaultContours (newdata, _sizex, _sizey, &fff);
    if (verbose)
      printf ("Contour 0B: Subsampling with %f\n", ContourSampling);
    //::SaveFits (newdata, _sizex, _sizey, "contour0b.fits");
    int nnx=0, nny=0;
    float *newdata1 =  LinearSampling (newdata, _sizex, _sizey, 
				       lineartrans, &nnx, &nny);
    float scale_x = (float)(_sizex)/(float)(nnx);
    float scale_y = (float)(_sizey)/(float)(nny);
    // Not sure why the xref is not correct for the contours...
//     float tr[6] = {_xref - _pixelscale_x, - _pixelscale_x*scale_x/1.25, 0.0, 
//     		   -_yref - _pixelscale_x/1.25, 
//     		   0.0, _pixelscale_y*scale_y,};     
    float tr[6] = {_xref- _pixelscale_x, - _pixelscale_x*scale_x, 0.0, 
		   -_yref + _pixelscale_y, 
		   0.0, _pixelscale_y*scale_y,};     

    //::SaveFits (newdata1, nnx, nny, "contour0c.fits");
    if (verbose)
      printf ("Contour 0C: Old size %d %d\n"
	      "            New size %d %d\n", _sizex, _sizey, nnx, nny);
    cpgconb (newdata1, nnx, nny, 1, nnx, 1, nny, 
	     &fff[0], -(int)fff.size(), tr, -1e30);
    delete[] newdata;
    delete[] newdata1;
  }
  else if (ContourSmoothing != 0.0) {
    double sigma = ContourSmoothing;
    if (verbose)
      printf ("Contour 1A: Smoothing with sigma=%.2f pixels\n", sigma);
    float *newdata =  FastSmooth (_data, _sizex, _sizey, sigma);
    float tr[6] = {_xref, -_pixelscale_x, 0.0, -_yref, 0.0, _pixelscale_y};
    if (fff.size()==0)
      DefaultContours (newdata, _sizex, _sizey, &fff);
    cpgcont (newdata, _sizex, _sizey, 1, _sizex, 1, _sizey, 
	     &fff[0], -(int)fff.size(), tr);
    delete[] newdata;
  }
  else if (ContourSampling != 1.0) {
    //float scale = ContourSampling;
    if (verbose)
      printf ("Contour 1B: Subsampling with %f\n", ContourSampling);
    int nnx=0, nny=0;
    if (fff.size()==0)
      DefaultContours (_data, _sizex, _sizey, &fff);
    float *newdata =  LinearSampling (_data, _sizex, _sizey, lineartrans, &nnx, &nny);
    if (verbose)
      printf ("%i %i -> %i %i, %f %f, %f \n", _sizex, _sizey, nnx, nny, _xref, _yref,_pixelscale_x);
    
    float scale_x = (float)(_sizex)/(float)(nnx);
    float scale_y = (float)(_sizey)/(float)(nny);
    // Not sure why the xref is not correct for the contours...
    float tr[6] = {_xref-_pixelscale_x, -_pixelscale_x*scale_x, 0.0, 
		   -_yref+_pixelscale_y, 
		   0.0, _pixelscale_y*scale_y,};     
    
    cpgconb (newdata, nnx, nny, 1, nnx, 1, nny, 
	     &fff[0], -(int)fff.size(), tr, -1e30);
    delete[] newdata;
  }
  else {
    float tr[6] = {_xref, -_pixelscale_x, 0.0, -_yref, 0.0, _pixelscale_y};
    if (verbose)
      printf ("Contour 1C: %i %i %f %f \n", _sizex, _sizey, _xref, _yref);
    if (fff.size()==0)
      DefaultContours (_data, _sizex, _sizey, &fff);
    //       cpgcont (data, sizex, sizey, 1, sizex, 1, sizey, 
    // 	       &fff[0], -(int)fff.size(), tr);
    cpgconb (_data, _sizex, _sizey, 1, _sizex, 1, _sizey, 
	     &fff[0], -(int)fff.size(), tr, -1e30);
  }



//   }
//   else {
//     int index = FindImage (lua_tostring (L,1));
//     if (index <0) {
//       const char *ch = lua_tostring (L,1);
//       std::cout << "ERROR: Cannot find image " << ch << std::endl;
//       lua_pushnumber (L, 1);
//       return (1);
//     }
//     if (ContourSmoothing != 0.0) {
//       double sigma = ContourSmoothing;
//       if (verbose)
// 	printf ("Smoothing with %f\n", sigma);
//       float *newdata =  FastSmooth (images[index]->data, 
// 				    images[index]->sizex, 
// 				    images[index]->sizey, sigma);
//       float scale = (float)sizex/(float)images[index]->sizex;
//       float tr[6] = {xref, -pixelscale_x*scale, 0.0, -yref, 
// 		     0.0, pixelscale_y*scale,};
//       if (fff.size()==0)
// 	DefaultContours (newdata, images[index]->sizex*images[index]->sizey, &fff);
//       cpgconb (newdata, images[index]->sizex, images[index]->sizey, 
// 	       1, images[index]->sizex, 1, images[index]->sizey, 
// 	       &fff[0], -(int)fff.size(), tr, -1e30);

//       delete[] newdata;
//     }
//     else if (ContourSampling != 1.0) {
//       float scale2 = (float)sizex/(float)images[index]->sizex;
//       if (verbose)
// 	printf ("Subsampling with %f\n", ContourSampling);
//       int nnx=0, nny=0;
//       float *newdata =  LinearSampling (images[index]->data, 
// 					images[index]->sizex, images[index]->sizey, 
// 					lineartrans, &nnx, &nny);
//       if (verbose)
// 	printf ("%i %i -> %i %i, %f %f \n", images[index]->sizex, images[index]->sizey, 
// 		nnx, nny, xref, yref);
      
//       float scale_x = (float)(images[index]->sizex)/(float)(nnx)*scale2;
//       float scale_y = (float)(images[index]->sizey)/(float)(nny)*scale2;
//       float tr[6] = {xref, -pixelscale_x*scale_x, 0.0, -yref, 
// 		     0.0, pixelscale_y*scale_y,};     
//       printf ("%f %f %f\n", scale2, scale_x, scale_y);
      
//       if (fff.size()==0)
// 	DefaultContours (newdata, nnx*nny, &fff);
//       cpgconb (newdata, nnx, nny, 1, nnx, 1, nny, 
// 	       &fff[0], -(int)fff.size(), tr, -1e30);
//       delete[] newdata;
//     }
//     else {
//       float scale = (float)sizex/(float)images[index]->sizex;
//       float tr[6] = {xref, -pixelscale_x*scale, 0.0, -yref, 
// 		     0.0, pixelscale_y*scale,};
//       if (fff.size()==0)
// 	DefaultContours (images[index]->data, images[index]->sizex*images[index]->sizey, &fff);
//       cpgconb (images[index]->data, images[index]->sizex, images[index]->sizey, 
// 	       1, images[index]->sizex, 1, images[index]->sizey, 
// 	       &fff[0], -(int)fff.size(), tr, -1e30);
//     }
// }
  if (lua_gettop(L) == 5)
    delete[] _data;
  return (Success(L));
}

// int Lua::SubSample (lua_State *L)
// {
//   int pad = 6, sfpad;
//   double lineartrans[6];
//   double scale;

//   if (lua_gettop(L) != 3)
//     return (MethodError (L, "Contour"));
  
//   int index = FindImage (lua_tostring (L,1));
//   if (index <0) {
//     const char *ch = lua_tostring (L,1);
//     std::cout << "ERROR: Cannot find image " << ch << std::endl;
//     lua_pushnumber (L, 1);
//     return (1);
//   }

//   scale = lua_tonumber (L,2);
//   sfpad = (int)(sf*pad);

//   in_ecl = image_new (images[index]->sizex+2*pad, images[index]->sizey+2*pad);
//   for (int i =0; i<images[index]->sizex; i++)
//     for (int j =0; j<images[index]->sizey; j++)
//       in_ecl->data[i+pad+(j+pad)*(sizex+2*pad)] 
// 	= images[index]->data[i+j*images[index]->sizex];
 
//   lineartrans[0] = scale ;
//   lineartrans[1] = 0.00 ;
//   lineartrans[2] = 0.00 ;
//   lineartrans[3] = 0.00 ;
//   lineartrans[4] = scale ;
//   lineartrans[5] = 0.00 ;
  
//   out_ecl = image_warp_linear(in_ecl, lineartrans, NULL);

//   Image *im = new Image;
//   if (limits) {
//     im->sminx = (int)(ch.val[0]+0.5);
//     im->sminy = (int)(ch.val[1]+0.5);
//     im->smaxx = (int)(ch.val[2]+0.5);
//     im->smaxy = (int)(ch.val[3]+0.5);
//   }
//   im->name = new char[strlen(ch.text[1])+1];
//   strcpy (im->name, ch.text[1]);

//   for (int i =sfpad; i<im_out->lx-sfpad; i++)
//     for (int j =sfpad; j<im_out->ly-sf*pad; j++)
//       temp[(i-sfpad)+(j-sfpad)*(im_out->lx-2*sfpad)] 
// 	= im_out->data[i+j*im_out->lx];


//   image_del (in_ecl);
//   image_del (out_ecl);
// }

//   static int SmoothGaussian (lua_State *L);
//int Lua::ClipImage (lua_State *L) 
//{
//  if (lua_gettop(L) != 6 || !lua_isstring(L,1) || !lua_isnumber(L,2) 
//      || !lua_isnumber(L,3) || !lua_isnumber(L,4) || !lua_isnumber(L,5) || !lua_isstring(L,6) )
//    return (MethodError (L, "Contour"));
//  Image *im = new Image;
//
//
//  im->sminx = (int)(ch.val[0]+0.5);
//  im->sminy = (int)(ch.val[1]+0.5);
//  im->smaxx = (int)(ch.val[2]+0.5);
//  im->smaxy = (int)(ch.val[3]+0.5);
//
//  im->name = new char[strlen(ch.text[1])+1];
//  strcpy (im->name, ch.text[1]);
// 
//}


int Lua::PxArith (lua_State *L)
{
  bool im;
  float val;
  if (lua_gettop(L) != 4) {
    std::cout << "ERROR: please specify exactly 4 arguments" 
	      << std::endl;
    return (MethodError (L, "ImArith"));
  }

  if (!(lua_isstring (L, 1) && lua_isstring(L,2) && lua_isstring (L,4)))
    return (MethodError (L, "ImArith"));

  int index1 = FindImage (lua_tostring(L,1));
  if (index1 <0) {
    std::cout << "ERROR: Cannot find image " << lua_tostring(L,1)
	      << std::endl;
    lua_pushnumber (L, 1);
    return (1);
  }
  //int index2 = FindImage (lua_tostring(L,4));

  if (lua_isnumber (L, 3)) {
    im = false;
    val = lua_tonumber (L, 3);
  }
  else{
  }

  printf ("THIS IS DUMMY\n");
  return (Success(L));
}

int Lua::PxReplace (lua_State *L)
{
  CheckArgsAlphaNumeric ch(L, 1, 3);
  if (!ch.Success())
    return (MethodError (L, "PxReplace"));

  int index = FindImage (ch.text[0]);
  if (index <0) {
    std::cout << "ERROR: Cannot find image " << ch.text[0] 
	      << std::endl;
    lua_pushnumber (L, 1);
    return (1);
  }

  for (unsigned int i=0; i<images[index]->sizex; i++)
    for (unsigned int j=0; j<images[index]->sizey; j++) {
      float val = images[index]->data[i+j*images[index]->sizex];
      if (val> ch.val[0] && val<ch.val[1])
	images[index]->data[i+j*images[index]->sizex] = ch.val[2];
    }
  return (Success(L));
}

int Lua::FreeImage (lua_State *L)
{
  CheckArgsAlphaNumeric ch(L, 1, 0);
  if (!ch.Success())
    return (MethodError (L, "FreeImage"));

  int index = FindImage (ch.text[0]);
  if (index <0) {
    std::cout << "ERROR: Cannot find image " << ch.text[0] 
	      << std::endl;
    lua_pushnumber (L, 1);
    return (1);
  }

  delete images[index];
  images.erase (images.begin()+index);
  return (Success(L));
}

int Lua::LoadImage (lua_State *L)
{
  int status = 0; 
  fitsfile *fptr;
  int hdutype;
  bool success;
  bool limits = false;
  int nnum;

  if (lua_gettop(L) == 2) 
    nnum = 0;
  else if (lua_gettop(L) == 6){
    nnum = 4;
    limits = true;
  }
  else
    return (MethodError (L, "LoadImage"));
  
  CheckArgsAlphaNumeric ch(L, 2, nnum);
  if (!ch.Success())
    return (MethodError (L, "LoadImage"));

  Image *im = new Image;
  if (limits) {
    im->sminx = (int)(ch.val[0]+0.5);
    im->sminy = (int)(ch.val[1]+0.5);
    im->smaxx = (int)(ch.val[2]+0.5);
    im->smaxy = (int)(ch.val[3]+0.5);
  }
  im->name = new char[strlen(ch.text[1])+1];
  strcpy (im->name, ch.text[1]);

  if ( !fits_open_file(&fptr, ch.text[0], READONLY, &status) ) {
    if (fits_get_hdu_type(fptr, &hdutype, &status) || hdutype != IMAGE_HDU) { 
      fits_close_file (fptr, &status);
      success = LoadDatImage (ch.text[0], im);
    }
    
    else {
      success = LoadFitsImage (fptr, im);
      fits_close_file (fptr, &status);
    }
  }
  else 
    success = LoadDatImage (ch.text[0], im);

  if (success) {
    images.push_back (im);
    return (Success (L));
  }

  delete im;
  std::cout << "ERROR: cannot open file " << ch.text[0] << std::endl;
  lua_pushnumber (L, 1);
  
  return (1);

}




int Lua::LoadFitsImage(fitsfile*fp, Image*image)
{
  int xmin=65536, xmax=-65536,ymin=65536, ymax=-65536;
  float *pix;
  int naxis, status = 0;
  long naxes[2],fpixel[2];

  fits_get_img_dim(fp, &naxis, &status);
  fits_get_img_size(fp, 2, naxes, &status);
  
  if (status || naxis != 2) { 
    printf("Error: NAXIS = %d.  Only 2-D images are supported.\n", naxis);
    abort ();
  }

  if (image->sminx == -65536)
    image->sminx = 1;
  if (image->sminy == -65536)
    image->sminy = 1;
  if (image->smaxx == -65536)
    image->smaxx = naxes[0];
  if (image->smaxy == -65536)
    image->smaxy = naxes[1];

  // Sanity checks
  if (image->sminx < 1 || image->sminx>naxes[0])
    image->sminx = 1;
  if (image->sminy < 1 || image->sminy>naxes[1])
    image->sminy = 1;
  if (image->smaxx < 1 || image->smaxx>naxes[0])
    image->smaxx = naxes[0];
  if (image->smaxy < 1 || image->smaxy>naxes[1])
    image->smaxy = naxes[1];
  if (image->sminx>=image->smaxx || image->sminy>=image->smaxy ) {
    printf ("Illegal figure x,y limits: using defaults\n");
    image->sminx = 1;
    image->sminy = 1;
    image->smaxx = naxes[0];
    image->smaxy = naxes[1];
  }

  image->sizex = image->smaxx-image->sminx+1;
  image->sizey = image->smaxy-image->sminy+1;
  image->xref = image->xref-image->sminx+1;
  image->yref = image->yref-image->sminy+1;

  image->data = new float[image->sizex*image->sizey+10000];
  fpixel[0] = image->sminx;
  pix = new float[naxes[0]];

  image->flmax = -1e10; image->flmin=1e10;
  for (fpixel[1] = image->sminy; fpixel[1] <= image->smaxy; fpixel[1]++) {
    if (fits_read_pix (fp, TFLOAT, fpixel, image->smaxx, 0, pix, 0, &status))
      break;;
    
    for (unsigned int i = 0; i<image->sizex; i++) {
      if (image->flmin>pix[i]) image->flmin=pix[i];
      if (image->flmax<pix[i]) image->flmax=pix[i];
      image->data[ i + image->sizex*(fpixel[1]-image->sminy+1)] = pix[i];
    }
  }
  delete[] pix;

  xmin = 1;
  xmax = image->sizex;
  ymin = 1;
  ymax = image->sizey;
  
  if (image->xref1 != 0.0 && image->yref1 != 0.0) {
    image->xref = (float)(image->xref1-(float)xmin)*image->pixelscale;
    image->yref = (float)(image->yref1-(float)ymin)*image->pixelscale;
  }
  else {
    image->xref = 0.5*(float)image->sizex*image->pixelscale;
    image->yref = 0.5*(float)image->sizey*image->pixelscale;
  }
  fits_close_file (fp, &status);
  return (1);
}

int Lua::LoadDatImage(char const*c, Image* image)
{
  std::ifstream is(c);

  std::vector<int> x;
  std::vector<int> y;
  std::vector<float> f;
  int xmin=65536, xmax=-65536,ymin=65536, ymax=-65536;

  if (is == NULL) {
    fprintf (stderr, "Cannot open file %s\n",c);
    exit (-1);
  }

  float  tmp1f, tmp2f, tmp3;

  while ( (is >> tmp1f) && (is >> tmp2f) && (is >> tmp3)  ) {
    x.push_back ((int)(tmp1f+0.1));
    y.push_back ((int)(tmp2f+0.1));
    //if (reverse) 
    // f.push_back (-tmp3);
    //else
    f.push_back (tmp3);
  }
  is.close ();

  for (unsigned int i = 0; i<x.size(); i++) {
    if (xmin>x[i]) xmin = x[i];
    if (xmax<x[i]) xmax = x[i];
    if (ymin>y[i]) ymin = y[i];
    if (ymax<y[i]) ymax = y[i];
    if (image->flmin>f[i]) image->flmin = f[i];
    if (image->flmax<f[i]) image->flmax = f[i];
  }

  if (image->smaxx > -65536 && image->smaxy > -65536 && 
      image->sminx > -65536 && image->sminy > -65536) {
    image->sizex = image->smaxx-image->sminx+1;
    image->sizey = image->smaxy-image->sminy+1;
    xmin = image->sminx;
    ymin = image->sminy;
    xmax = image->smaxx;
    ymax = image->smaxy;
  }
  else if (image->smaxx > -65536 && image->smaxy > -65536) {
    image->sizex = image->smaxx;
    image->sizey = image->smaxy;
    xmax = image->smaxx;
    ymax = image->smaxy;
    xmin = 1;
    ymin = 1;
  }
  else {
    image->sizex = xmax - xmin+1;
    image->sizey = ymax - ymin+1;
  }

  image->data = new float[ image->sizex*image->sizey ];    

  for (unsigned int i = 0; i<image->sizex*image->sizey; i++)
    image->data[i] = -1e30;

  for (unsigned int i = 0; i<x.size(); i++) {
    if (x[i]>=xmin & x[i]<=xmax & y[i]>=ymin & y[i]<=ymax)
      image->data[ (x[i]-xmin) + image->sizex*(y[i]-ymin)] = f[i];
  }

  if (image->xref1 != 0.0 && image->yref1 != 0.0) {
    image->xref = (float)(image->xref1-(float)xmin)*image->pixelscale;
    image->yref = (float)(image->yref1-(float)ymin)*image->pixelscale;
  }
  else {
    image->xref = 0.5*(float)image->sizex*image->pixelscale;
    image->yref = 0.5*(float)image->sizex*image->pixelscale;
  }
  return (1);
}

int Lua::MethodError (lua_State *L, const char *method)
{
  std::cout << "ERROR: incorrect input in " << method << std::endl;
  lua_pushnumber (L, 1);
  return (1);
}

int Lua::Success (lua_State *L)
{
  lua_pushnumber (L, 0);
  return (1);
}

int Lua::ViewPort (lua_State *L)
{
  CheckArgsNumericAlpha ch(L, 4, 0);
  if (!ch.Success())
    return (MethodError (L, "ViewPort (x1, y1, x2, y2)"));

  cpgsvp (ch.val[0],ch.val[2], ch.val[1], ch.val[3]);

  lua_pushnumber (L, 0);
  return (1);
}

int Lua::Limits (lua_State *L)
{
  CheckArgsNumericAlpha ch(L, 4, 0);
  if (!ch.Success())
    return (MethodError (L, "Limits"));

  cpgswin (ch.val[0],ch.val[1], ch.val[2], ch.val[3]);

  lua_pushnumber (L, 0); 
  return (1);
}

int Lua::Plot (lua_State *L)
{
  if (lua_gettop(L)!=4) {
    printf ("Usage: \nPlot (\"file\", column1, column2, \"type\");\n");
    return (MethodError (L, "Plot"));
  }

  if (! (lua_isstring(L, 1) && lua_isnumber(L, 2) &&
	 lua_isnumber(L, 3) && lua_isstring(L, 4) )) {
    printf ("Usage: \nPlot (\"file\", column1, column2, \"type\");\n");
    return (MethodError (L, "Plot"));
  }

  const char *fname = lua_tostring (L, 1);
  int c1 = (int)(lua_tonumber (L, 2)+0.5);
  int c2 = (int)(lua_tonumber (L, 3)+0.5);
  const char *type = lua_tostring (L, 4);

  struct stat st;
  if (stat (fname, &st) == -1) {
    perror("Plot: stat");
    return (MethodError (L, "Plot"));
  }
  FILE *file;
  file = fopen (fname, "r");
  char *buf = new char[st.st_size];
  fread (buf, st.st_size, 1, file);
  fclose (file);

  std::vector<float> x;
  std::vector<float> y;

  int pos = 0;
  while (pos<st.st_size) {
    int pos1 = pos;
    while (buf[pos1]!='\n')
      pos1++;

    int i = 0;
    
    while (buf[pos]==' ' || buf[pos]=='\t')
      pos++;

    float tx = atof (&buf[pos]);
    while (buf[pos]!=' ' && buf[pos]!='\t')
      pos++;
    while (buf[pos]==' ' || buf[pos]=='\t')
      pos++;
    float ty = atof (&buf[pos]);
    while (buf[pos]!='\n')
      pos++;
    if (pos != pos1)
      pos = pos1+1;
    else {
      pos = pos1 + 1;
      x.push_back (tx);
      y.push_back (ty);
    }
  }

  if (x.size()>0)
    cpgline (x.size(), &x[0], &y[0]);
  delete[] buf;

  lua_pushnumber (L, 0); 
  return (1);
}

int Lua::Box (lua_State *L)
{
  if (lua_gettop(L)!=6) {
    printf ("Usage: \nBox (xopt, xtick, nxsub, yopt, ytick, nysub)\n");
    return (MethodError (L, "Limits"));
  }
  if (! (lua_isstring(L, 1) && lua_isnumber(L, 2) &&
	 lua_isnumber(L, 3) && lua_isstring(L,4) &&
	 lua_isnumber(L, 5) && lua_isnumber (L, 6) ) ) {
    printf ("Usage: \nBox (xopt, xtick, nxsub, yopt, ytick, nysub)\n");
    return (MethodError (L, "Limits"));
  }
  
  const char *xopt = lua_tostring (L, 1);
  float xtick = lua_tonumber (L, 2);
  int nxsub = (int)(0.5+lua_tonumber(L, 3));
  const char *yopt = lua_tostring (L, 4);
  float ytick = lua_tonumber (L, 5);
  int nysub = (int)(0.5+lua_tonumber(L, 6));

  float x1, x2, y1, y2;
  int ci;
  cpgqwin (&x1, &x2, &y1, &y2);
  cpgqci (&ci);
  cpgsci (0);
  cpgrect (x1, x2, y1, y2);
  cpgsci (ci);
  cpgbox (xopt, xtick, nxsub, yopt, ytick, nysub);

  lua_pushnumber (L, 0);
  return (1);
}

int Lua::FindImage (const char *c)
{
  for (unsigned int i = 0; i<images.size(); i++) 
    if (strcmp(images[i]->name, c)==0)
      return (i);

  return (-1);
}
/*
 *   Classes used to 
 *     check argument types
 *     place the values to C-type array
 *   used internally in Lua 
 */
Lua::CheckArgsNumeric::CheckArgsNumeric (lua_State *L,int n)
  :val(NULL), success(true)
{
  if (n != lua_gettop(L)) {
    std::cout << "ERROR: please specify exactly " << n << " arguments" 
	      << std::endl;
    success = false;
  }
  else {
    int i=1;
    while (i<=n) {
      if (!lua_isnumber(L, i)) {
	std::cout << "ERROR: argument number " << i << " is not a number" 
		  << std::endl;
	success = false;
      }
      i++;
    }
    if (success) {
      val = new float[n];
      for (int i=0; i<n; i++)
	val[i] = lua_tonumber(L, i+1);
    }
  }
}

Lua::CheckArgsNumeric::~CheckArgsNumeric ()
{
  if (success)
    delete[] val;
}

Lua::CheckArgsNumericAlpha::CheckArgsNumericAlpha (lua_State *L,int nf, int ns)
  :val(NULL), text(NULL),success(true)
{
  if ((ns+nf) != lua_gettop(L)) {
    std::cout << "ERROR: please specify exactly " << (ns+nf) << " arguments" 
	      << std::endl;
    success = false;
  }
  else {
    int i=1;
    while (i<=nf) {
      if (!lua_isnumber(L, i)) {
	std::cout << "ERROR: argument number " << i << " is not a number" 
		  << std::endl;
	success = false;
      }
      i++;
    }
    i = 1;
    while (i<=ns) {
      if (!lua_isstring(L, i+nf)) {
	std::cout << "ERROR: argument number " << (nf+i) 
		  << " is not a string" << std::endl;
	success = false;
      }
      i++;
    }
    if (success) {
      val = new float[nf];

      for (int i=0; i<nf; i++)
	val[i] = lua_tonumber(L, i+1);

      text = new const char*[ns];
      for (int i=0; i<ns; i++)
	text[i] = lua_tostring(L, nf+i+1);
    }
  }
}

Lua::CheckArgsNumericAlpha::~CheckArgsNumericAlpha ()
{
  if (success) {
    delete[] text;
    delete[] val;
  }
}


Lua::CheckArgsAlphaNumeric::CheckArgsAlphaNumeric (lua_State *L,int ns, int nf)
  :val(NULL), text(NULL),success(true)
{
  if ((ns+nf) != lua_gettop(L)) {
    std::cout << "ERROR: please specify exactly " << (ns+nf) << " arguments" 
	      << std::endl;
    success = false;
  }
  else {
    int i=1;
    while (i<=ns) {
      if (!lua_isstring(L, i)) {
	std::cout << "ERROR: argument number " << i << " is not a string" 
		  << std::endl;
	success = false;
      }
      i++;
    }
    i = 1;
    while (i<=nf) {
      if (!lua_isstring(L, i+ns)) {
	std::cout << "ERROR: argument number " << (ns+i) 
		  << " is not a number" << std::endl;
	success = false;
      }
      i++;
    }
    if (success) {
      val = new float[nf];

      for (int i=0; i<nf; i++)
	val[i] = lua_tonumber(L, ns+i+1);

      text = new const char*[ns];
      for (int i=0; i<ns; i++)
	text[i] = lua_tostring(L, i+1);
    }
  }
}

Lua::CheckArgsAlphaNumeric::~CheckArgsAlphaNumeric ()
{
  if (success) {
    delete[] text;
    delete[] val;
  }
}

Image::Image ()
  : data(NULL), sizex(0), sizey(0), xref(-60000.0), yref(-60000.0),
    pixelscale (1.0), name(NULL),
    sminx(-65536), smaxx(-65536), sminy(-65536), smaxy(-65536),
    xref1(0.0),yref1(0.0)
{

}


Image::~Image ()
{
  if (data != NULL)
    delete[] data;
  if (name != NULL)
    delete[] name;
}



/*
 *  PGPLOT HELPER ROUTINES!
 */ 



void pgellipse_c( double Cpx, 
                  double Cpy, 
                  double Major, 
                  double Minor, 
                  double Pa,
                  double Startangle,
                  double Endangle,
                  double Delta )
/*------------------------------------------------------------*/
/* Draw ellipse from 'Startangle' to 'Endangle' */
/*------------------------------------------------------------*/
{
   double    CosP, SinP;                    /* Angles */
   double    CosA, SinA;                    /* Angles */
   double    Alpha;                         /* Used in Polar coordinates */
   double    R;                             /* Radius used in Polar coordinates */
   double    Denom;                         /* Help var. */
   double    Xell, Yell;                    /* Points of not rotated ellipse */
   //double    Xrot, Yrot;                    /* X,Yell rotated over Pos.angle */
   float     *Xpoints, *Ypoints;            /* Outline ellipse in arrays */
   int       i;                             /* Array index */
   int       len;                           /* Array length */
   int       Numpoints;                     /* Number of points in arrays */
   
   CosP  = cos( Pa * 0.01745329251994329576923 );
   SinP  = sin( Pa * 0.01745329251994329576923 );   

   len = (int) ( fabs((Endangle)-(Startangle))/(Delta) );   
   len += 10;                                            /* to be (very) sure */
   Xpoints = (float *) calloc( len, sizeof(float) );
   Ypoints = (float *) calloc( len, sizeof(float) );   

   
   if (Xpoints == NULL || Ypoints == NULL)
   {
      return;
   }


   i = 0;
   for (Alpha  = (double) Startangle; 
        Alpha <= (double) Endangle; 
        Alpha += (double) Delta) 
   {
      /*-----------------------------------------------------------*/
      /* Ellipse: b^2.x^2 + a^2.y^2 = a^2.b^2                      */
      /* Substitute: x = r.cos(alpha), y = r.sin(alpha) and solve  */
      /* for r. Repeat this action for angles between 0 and 90 deg */
      /*-----------------------------------------------------------*/      
      CosA = cos(Alpha* 0.017453292519943295769237);
      SinA = sin(Alpha* 0.017453292519943295769237);      
      Denom = (Minor*CosA * Minor*CosA + Major*SinA * Major*SinA);
      if (Denom == 0.0) {
         R = 0;
      }
      else {
         R = sqrt( Minor*Major * Minor*Major / Denom );
      }   
      Xell = R * CosA;
      Yell = R * SinA;
      /* We have a point on this ellipse, now rotate this point */
      /* and move to given origin */
      Xpoints[i] = (float) (Xell * CosP - Yell * SinP  + Cpx); 
      Ypoints[i] = (float) (Xell * SinP + Yell * CosP  + Cpy);
      i++;
   }
   
   Numpoints = i;
   cpgline( Numpoints, Xpoints, Ypoints );   /* Can be replaced by 'pgpoly' */         
   free( Ypoints );
   free( Xpoints );
      
} 

void pgcone_c( double Cpx, 
	       double Cpy, 
	       double R, 
	       double Startangle,
	       double Endangle,
	       double Delta )
{
   double    CosA, SinA;                    /* Angles */
   double    Alpha;                         /* Used in Polar coordinates */
   float     *Xpoints, *Ypoints;            /* Outline ellipse in arrays */
   int       i;                             /* Array index */
   int       len;                           /* Array length */
   int       Numpoints;                     /* Number of points in arrays */
   

   len = (int) ( fabs((Endangle)-(Startangle))/(Delta) );   
   len += 10;                                            /* to be (very) sure */
   Xpoints = (float *) calloc( len, sizeof(float) );
   Ypoints = (float *) calloc( len, sizeof(float) );   

   
   if (Xpoints == NULL || Ypoints == NULL)
   {
      return;
   }


   Xpoints[0] = Cpx;
   Ypoints[0] = Cpy;

   i = 1;
   for (Alpha  = (double) Startangle; 
        Alpha <= (double) Endangle; 
        Alpha += (double) Delta) 
   {
      CosA = cos(Alpha* 0.017453292519943295769237);
      SinA = sin(Alpha* 0.017453292519943295769237);      

      Xpoints[i] = (float) (Cpx + R*SinA); 
      Ypoints[i] = (float) (Cpy + R*CosA);
      i++;
   }
   
   Xpoints[i] = Cpx;
   Ypoints[i] = Cpy;
   i++;

   Numpoints = i;
   cpgline( Numpoints, Xpoints, Ypoints );   /* Can be replaced by 'pgpoly' */         
   free( Ypoints );
   free( Xpoints );
      
} 



float *  LinearSampling (float *array, int nx, int ny, 
			 double *trans, int *nnx, int *nny)
{
  image_t * im_in, * im_out;
  float *temp;
  float sf = 0.5*(trans[0]+trans[4]);

  int pad = 0;
  int sfpad = (int)(sf*pad);
  int isf = (int)(sf+0.1);
  if (verbose)
    printf ("LinearSampling:\n  pad: %d\n sfpad: %d\n isf: %d\n", 
	    pad, sfpad, isf);
  
  im_in = image_new (nx+2*pad, ny+2*pad);

  for (int i =0; i<nx; i++)
    for (int j =0; j<ny; j++)
      im_in->data[i+pad+(j+pad)*(nx+2*pad)] = array[i+j*nx];
  
  for (int i =0; i<nx; i++)
    for (int j =0; j<ny; j++)
      if (array[i+j*nx]<-1e29) {
	float ff = 0.0;
	float nn = 0.0;
	if (i-1>=0 && array[i-1+j*nx]>-1e29) {
	  ff = ff + array[i-1+j*nx];
	  nn = nn + 1.0;
	}
	if (i+1<nx && array[i+1+j*nx]>-1e29) {
	  ff = ff + array[i+1+j*nx];
	  nn = nn + 1.0;
	}

	if (i-1>=0 && j-1>=0 && array[i-1+(j-1)*nx]>-1e29) {
	  ff = ff + array[i-1+(j-1)*nx];
	  nn = nn + 1.0;
	}
	if (j-1>=0 && array[i+(j-1)*nx]>-1e29) {
	  ff = ff + array[i+(j-1)*nx];
	  nn = nn + 1.0;
	}
	if (j+1<nx && j-1>=0 && array[i+1+(j-1)*nx]>-1e29) {
	  ff = ff + array[i+1+(j-1)*nx];
	  nn = nn + 1.0;
	}

	if (i-1>=0 && j+1<ny && array[i-1+(j+1)*nx]>-1e29) {
	  ff = ff + array[i-1+(j+1)*nx];
	  nn = nn + 1.0;
	}
	if (j+1<ny && array[i+(j+1)*nx]>-1e29) {
	  ff = ff + array[i+(j+1)*nx];
	  nn = nn + 1.0;
	}
	if (j+1<nx && j+1<ny && array[i+1+(j+1)*nx]>-1e29) {
	  ff = ff + array[i+1+(j+1)*nx];
	  nn = nn + 1.0;
	}
	if (nn>0.5)
	  im_in->data[i+pad+(j+pad)*(nx+2*pad)] = ff/nn;
	else
	  im_in->data[i+pad+(j+pad)*(nx+2*pad)] = 0.0;
      }
      else
	im_in->data[i+pad+(j+pad)*(nx+2*pad)] = array[i+j*nx];

//   for (int k=0; k<pad; k++) {
//     for (int i=pad;i<nx+pad; i++) {
//       im_in->data[i+k*(nx+2*pad)]=im_in->data[i+pad*(nx+2*pad)];
//       im_in->data[i+(ny+2*pad-1-k)*(nx+2*pad)]=im_in->data[i+(ny+pad-1)*(nx+2*pad)];
//     }
    
//     for (int j=pad;j<ny+pad; j++) {
//       im_in->data[k+j*(nx+2*pad)] = im_in->data[pad+j*(nx+2*pad)];
//       im_in->data[nx+2*pad-1-k+j*(nx+2*pad)] = im_in->data[nx+pad-1+j*(nx+2*pad)];
//     }
//   }

  if (verbose)
    printf ("Warping\n");
  im_out = image_warp_linear(im_in, trans, NULL);
  if (verbose)
    printf ("Warping done\n");

  temp = new float[(im_out->lx-2*sfpad)*(im_out->ly-2*sfpad)];

  for (int i =sfpad; i<im_out->lx-sfpad; i++)
    for (int j =sfpad; j<im_out->ly-sf*pad; j++)
      temp[(i-sfpad)+(j-sfpad)*(im_out->lx-2*sfpad)] 
	= im_out->data[i+j*im_out->lx];

  for (int i =0; i<nx; i++)
    for (int j =0; j<ny; j++) 
      if (array[i+j*nx]<-1e29) {
	for (int ii = 0; ii<isf; ii++)
	  for (int jj = 0; jj<isf; jj++) {
	    temp[ii+i*isf+(jj+j*isf)*(im_out->lx-2*sfpad)]
	      = -1e30;
	  }
      }

//   for (int i =0; i<nx; i++)
//     for (int j =0; j<ny; j++) {
//       int nei=0;
//       for (int ii=-1; ii<2; ii++)
// 	for (int jj=-1; jj<2; jj++)
// 	  if ( (i+ii)>=0 && (i+ii)<nx
// 	       && (j+jj)>=0 && (j+jj)<ny
// 	       && (jj!=0 && ii!=0)) 
// 	    if (array[i+ii+(j+jj)*nx]>-1e29)
// 	      nei++;
//       if (nei < 2) {
// 	for (int ii = 0; ii<isf; ii++)
// 	  for (int jj = 0; jj<isf; jj++) 
// 	    temp[ii+i*isf+(jj+j*isf)*(im_out->lx-2*sfpad)]
// 	      = array[i+j*nx];
	
//       }
//     }
    
  *nnx = im_out->lx-2*sfpad;
  *nny = im_out->ly-2*sfpad;

  image_del (im_in);
  image_del (im_out);

  return (temp);
}

float *Smooth (float *array, int nx, int ny, double sigma)
{
  int halfbox = (int)(4.0*sigma+1.0);
  int fullbox = 2*halfbox+1;
  float *arraysum = new float[nx*ny];
  float *array2 = new float[nx*ny];

  float *kernel = new float[fullbox*fullbox];
  float sum = 0.0;

  if (verbose)
    printf ("halfbox = %d\n",halfbox);

  for (int i=-halfbox; i<=halfbox; i++)
    for (int j=-halfbox; j<=halfbox; j++) {
      float rad = sqrt(pow(i, 2.0)+pow(j, 2.0));
      kernel[i+halfbox+(j+halfbox)*fullbox]
	= exp(-0.5*pow(rad/sigma, 2.0));
      sum = sum + kernel[i+halfbox+(j+halfbox)*fullbox];
    }
  
  for (int i=0; i<fullbox; i++)
    for (int j=0; j<fullbox; j++)
      kernel[i+j*fullbox]=kernel[i+j*fullbox]/sum;

  
  for (int i = 0; i<nx; i++)
    for (int j = 0; j<ny; j++) {
      arraysum[i+j*nx] = 0.0;
      array2[i+j*nx] = 0.0;
    }

  for (int i = 0; i<nx; i++)
    for (int j = 0; j<ny; j++) 
      for (int ii=-halfbox; ii<=halfbox; ii++)
	for (int jj=-halfbox; jj<=halfbox; jj++) 
	  if ( (i+ii)>0 && (j+jj)>0 && (i+ii)<nx && (j+jj)<ny) {
	    array2[i+ii+(j+jj)*nx] = array2[i+ii+(j+jj)*nx] +
	      kernel[ii+halfbox+(jj+halfbox)*fullbox]
	      * array[i+j*nx];
	      arraysum[i+ii+(j+jj)*nx] = arraysum[i+ii+(j+jj)*nx]
		+ kernel[ii+halfbox+(jj+halfbox)*fullbox];
	  }

  // Final normalization (in practise, this will only affect the edges)
  for (int i = 0; i<nx; i++)
    for (int j = 0; j<ny; j++) 
      array2[i+j*nx] = array2[i+j*nx] / arraysum[i+j*nx];

  delete[] kernel;
  delete[] arraysum;

  // SaveFits (array2, nx, ny, "testislow.fits");
  return (array2);

}

void SaveFits (float *array, int nx, int ny, const char*name)
{
  fitsfile *fptr;
  int status = 0;

  long naxes[2] = {nx, ny};
  long fpixel[2] = {1, 1};
  fits_create_file  (&fptr, name , &status);
  if (status) {
    std::cout << "ERROR: Cannot create file " << name
	      << " - it probably exists already" 
	      << std::endl;
    return;
  }

  fits_create_img(fptr, FLOAT_IMG, 2, naxes, &status);

  fits_write_pix(fptr, TFLOAT, fpixel, nx*ny,array, &status);
  fits_close_file (fptr, &status);
}

float *FastSmooth (float *array, int nx, int ny, double sigma)
{
  int halfbox = (int)(4.0*sigma+1.0);
  int fullbox = 2*halfbox+1;
  float *arraysum = new float[nx*ny];
  float *array2 = new float[nx*ny];
  float *array3 = new float[nx*ny];

  float *kernel = new float[fullbox];
  float sum = 0.0;

  if (verbose)
    printf ("halfbox = %d\n",halfbox);

  for (int i=-halfbox; i<=halfbox; i++) {
    float rad = (float)i;
    kernel[i+halfbox] = exp(-0.5*pow(rad/sigma, 2.0));
    sum = sum + kernel[i+halfbox];
  }

  for (int i=0; i<fullbox; i++) {
    kernel[i]=kernel[i]/sum;
  }
    
  for (int i = 0; i<nx*ny; i++) {
    arraysum[i] = 0.0;
    array2[i] = 0.0;
    array3[i] = 0.0;
  }

  // First direction
  for (int j = 0; j<ny; j++) 
    for (int i = 0; i<nx; i++) {
      int ii1 = -halfbox;
      int ii2 = halfbox;
      if (i-halfbox <0) ii1 = -i;
      if (i+halfbox>=nx) ii2=nx-i-1;
      for (int ii=ii1; ii<=ii2; ii++) {
	array2[i+ii+(j)*nx] = array2[i+ii+(j)*nx] +
	  kernel[ii+halfbox]
	  * array[i+j*nx];
	arraysum[i+ii+(j)*nx] = arraysum[i+ii+(j)*nx]
	  + kernel[ii+halfbox];
      }
    }
     
  for (int i = 0; i<nx*ny; i++) {
    array2[i] = array2[i] / arraysum[i];
    arraysum[i] = 0.0;
  }

  // 2nd direction
  for (int i = 0; i<nx; i++)
    for (int j = 0; j<ny; j++) {
      int jj1 = -halfbox;
      int jj2 = halfbox;
      if (j-halfbox <0) jj1 = -j;
      if (j+halfbox>=ny) jj2=ny-j-1;
      for (int jj=jj1; jj<=jj2; jj++) {
	array3[i+(j+jj)*nx] = array3[i+(j+jj)*nx] +
	  kernel[jj+halfbox]
	  * array2[i+j*nx];
	arraysum[i+(j+jj)*nx] = arraysum[i+(j+jj)*nx]
	  + kernel[jj+halfbox];
      }
    }

  // Final normalization (in practise, this will only affect the edges)
  for (int i = 0; i<nx*ny; i++)
    array3[i] = array3[i] / arraysum[i];

  delete[] kernel;
  delete[] arraysum;
  delete[] array2;
  
  // SaveFits (array3, nx, ny, "testifast2.fits");
  return (array3);

}

void DefaultContours (float *array, int nx, int ny, std::vector<float> *vec)
{
  double sum = 0.0;
  double sum2 = 0.0;
  double n=0.0;
  double mean, rms, oldrms;
  int npix=nx*ny;
  for (int i=0; i<npix; i++)
    if (std::isfinite(array[i]) && array[i]!=0.0) {
      sum = sum+array[i];
      sum2= sum2+pow(array[i], 2.0);
      n = n+ 1.0;
    }
  mean = sum/n;
  rms = sqrt(sum2/n -pow(mean, 2.0));

  for (int k =1; k<10; k++) {
    int i;
    oldrms = rms;
    sum  = 0.0;
    sum2 = 0.0;
    n    = 0.0;
    for (int ii=0; ii<nx; ii++)
      for (int jj=0; jj<ny; jj++) {
	i = ii+nx*jj;
	if (std::isfinite(array[i]) && array[i]!=0.0 && fabs(array[i]-mean)<3.0*rms) {
	  sum  += array[i];
	  sum2 += array[i]*array[i];
	  n++;
	}
      }
    mean = sum/n;
    rms = sqrt(sum2/n -pow(mean, 2.0));
    if (rms>0.99*oldrms) {
      if (verbose)
	printf ("RMS estimation finished after %d iterations\n",k );
      break;
    }
  }
  
  printf ("DefaultContours: background estimation, mean: %f, rms: %f\n", mean, rms);

  vec->push_back (mean+ContourStart*rms);
  for (int k=0; k<32; k++)
    vec->push_back (mean +ContourStart*rms + 2.0*rms*pow(2.0,(float)k));

}


int Fits2Float (const char *imname, float **array, int *arx, int *ary)
{
  int naxis, status = 0;
  long naxes[3],fpixel[3];
  int hdutype, hdunum;
  fitsfile *fp;
  if (verbose)
    printf ("Fits2Float\n");

  naxes[0] = 1;
  naxes[1] = 1;
  naxes[2] = 1;
  
  if ( fits_open_file(&fp, imname, READONLY, &status) ) {
    printf ("Cannot open file %s\n", imname);
    fits_close_file (fp, &status);
    return (-1);
  }

  if (fits_get_hdu_type(fp, &hdutype, &status) || hdutype != IMAGE_HDU) { 
    printf ("Cannot open file %s: not an image\n", imname);
    fits_close_file (fp, &status);
    return (-1);
  }

  fits_get_img_dim(fp, &naxis, &status);
  fits_get_img_size(fp, naxis, &naxes[0], &status);

  if (naxis == 0) {
    fits_get_num_hdus (fp, &hdunum, &status);
    if (hdunum<2 || status) {
      printf("Error: HDUNUM = %d.\n", hdunum);
      abort ();
    }

    fits_movabs_hdu (fp, 2, &hdutype, &status);
    fits_get_img_dim(fp, &naxis, &status);
    fits_get_img_size(fp, naxis, &naxes[0], &status);
  }

  if (status || naxis != 2) { 
    printf("Error: NAXIS = %d.  Only 2D images are supported.\n", naxis);
    abort ();
  }
  
  if (status==1) {
    printf ("ERROR: CFitsio has returned an error in LoadFitsFile.\n"
	    "       Perhaps the file is corrupted?\n");
    abort();
  }
  
  *array = new float[naxes[0]*naxes[1]+10000];
  fpixel[0] = 1; fpixel[1] = 1; fpixel[2]= 1;
  
  fits_read_pix (fp, TFLOAT, fpixel, naxes[0]*naxes[1], 
		 0, (*array), 0, &status);
  if (status==1) {
    delete[] array;
    printf ("ERROR: fits_read_pix has returned an error in LoadFitsFile.\n"
	    "       Perhaps the file is corrupted?\n");
    abort();
  }
  (*arx) = naxes[0];
  (*ary) = naxes[1];
  fits_close_file (fp, &status);

  //printf ("Saving\n");
  //SaveFits (*array, naxes[0], naxes[1], "testislow.fits");
  return (0);
}
